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ABSTRACT 


An efficient method of computing structural response of multi-story nonlinear 
base isolated buildings for a given seismic event is presented. Using a recursive block- 
by-block integral equation formulation (RBBIEF) solution to the governing nonlinear 
Volterra integral equation, structural base motion acting through an arbitrary number of 
nonlinear base isolators can be computed quickly and accurately. The general solution to 
the governing nonlinear Volterra integral is formulated and subsequently converted into 
code using MATLAB. The procedure incorporates modal properties, computed from 
conventional finite element (FE) techniques, and the generated MATLAB programs to 
solve a varying set of multi-degree of freedom structures coupled to both linear and 
nonlinear isolators. Ultimately, an analysis is conducted on a 30-story building that was 
overly designed using the 1994 Load Resistance Factor Design [Ref. 1] and the 1994 
Uniform Building Codes [Ref. 2] for earthquake loading. The method demonstrates that 
the Volterra integration scheme in the time domain is very effective and efficient. 
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I. INTRODUCTION 


Seismic base motion is an important design factor in constructing and retrofitting 
structures. The extensive base displacements and resultant shear forces caused by large 
magnitude earthquakes are major contributors to the eventual failure of the structure. 
Isolating structures from seismic events is currently accomplished through the use of 
passive isolators. 

Originally, it was the engineer’s goal to reduce the structural displacements 
through isolation. However, current techniques place equal importance on increasing the 
energy dissipation in an isolation device. An increase in the area contained within the 
device’s hysteresis loop is directly related to the energy dissipation provided by the 
device. The plastic properties of passive hysteretic type isolators are mathematically 
modeled as a complex nonlinear equation [Ref. 3]. Modeling and analysis using finite 
element (FE) techniques for nonlinear isolators then becomes equally complex and 
computationally demanding. 

Several Commercial-Off-the-Shelf (COTS) programs such as, IDRAC-3D, 3D- 
BASIS, DRAIN-2DX, SAP, ANSR, N-PAD, LPM, and others have their respective 
limitations and advantages as outlined in [Ref. 3]. The primary solution methods for 
nonlinear base isolation of structures are the pseudo-force iterative method or the 
nonlinear tangent stiffness method [Ref. 3]. However, the incorporation of 3-D motion 
and discrete time histories in conjunction with n-nonlinear isolators to the above- 
mentioned FE programs can be computationally exhaustive and or not practically 
possible. 
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Current FE programs can, in seconds, compute normalized mode shapes and 
natural frequencies of several thousand degree of freedom systems. The structure’s 
modal properties are used by many of the above programs in order to compute transient 
response or frequency response. Unlike the methods of solving the differential equations 
of motion above, we describe herein a method that solves the transient response directly 
using the Volterra integral form. This method, referred to as time domain synthesis, is 
presented in [Ref. 4] and [Ref. 5], and is a numerically accurate and fast solution for 
stractural response due to base motion excitation through nonlinear isolation. 

The addition of nonlinear memory-type isolators to the time domain synthesis is 
presented herein. Currently, the synthesis is solved using the (RBBIEF) method. The 
RBBIEF method obtains its significant efficiency from two characteristics of the 
formulation. First, a stmcture is described by impulse response functions calculated only 
for those physical DOF of interest. Those DOF for which response is not required can be 
omitted with no loss of accuracy. Second, the block-by-block convolution divides a large 
time history into several “blocks,” and the governing Volterra integral equation is solved 
in a block-by-block manner wherein only one block undergoes iterative convolution at a 
time. The converged forces of synthesis from previously iterated blocks are used in 
conjunction with the results of the current iteration to calculate the response for the 
current block. To incorporate memory-type isolators in the synthesis, the isolator 
equations are solved using standard numerical solution techniques for ordinary 
differential equations. This solution of the isolator’s governing differential equation is 
incorporated into the RBBIEF method. This revised synthesis method provides fast and 



accurate solutions. Ultimately, the combination of quick eigenvalue solutions of large 
degree of freedom structures, modal reduction, and a block-by-block recursive solution 
of the governing integral equation of the retained nodes that are seismically excited, can 
lead to significant timesavings and accuracy. 
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11. EARTHQUAKE DESIGN AND NONLINEAR ISOLATORS 


A. STRUCTURAL DESIGN 

Current design techniques and survivability rely upon the structure’s total strain 
energy capacity as it plastically deforms during a seismic event. This “hard working” 
process is accounted for by the Uniform Building Code, UBC. In Figure 1 below, the 
UBC’s required seismic coefficient is significantly underestimating the various 
earthquakes presented. 



Figure 1. UBC Seismic Acceleration Graph From Ref. [3] 

Additional varying normalized response spectra plots similar to Figure 1 are 
produced in [Ref. 3]. The vertical axis outlines the seismic coefficient and is used in 
various UBC equations, the horizontal axis represents structural periods. The apparent 
differences between the UBC’s required protection as seen in Figure 1 and an arbitrary 
seismic excitation is where passive base isolation plays it’s important role. In order to 
ensure the reduction of the earthquake’s overall input forces using passive isolators, an 
assumed normalized response spectrum is chosen. 
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Existing codes outlined in Chapter 23 of [Ref. 3] require that the existing 
structure have several design capacities such as seismic base shear, lateral seismic force, 
lateral deflection, overturning moment, and torsional moment. As part of the seismic 
design procedure, the structure is converted into modal space and a conservative 
participation factor is chosen. The general design equation is presented as: 



Replacing \|/ with V for shear, A for drift, 5 for deflection, M for moment, and T for 
torsion into Equation (1) above to arrive at the respective design limits. The subscript 
(m) represents the modal values and (x) the story level. Details of modal decomposition 
and participation can be found in numerous vibration textbooks such as [Ref. 6] and will 
be summarized in subsequent chapters herein. 

A 30 story moment resisting frame (MRF) building was conservatively designed 
using the Manual of Steel Construction, Load Resistance Factor Design (LRFD) [Ref.l], 
the Uniform Building Code (UBC). [Ref.2], and a steel structure design textbook [Ref. 
7]. The structure has an inter-story height of 17 feet and a total building height of 510 ft. 
The columns were spaced 25-ft center to center leading to an overall length of 400 ft at 
the base. A total of 256 columns are present corresponding to 256 base isolators at the 
base. Each floor was designed to support a dead load of 150 psf and the roof, 50 psf. The 
beams were divided as follows: the roof consisted of W36xl35 beams; Floors 30-21 
W36xl50; Floors 20-11 W36xl60; Floors 10-2 W36xl70. The first floor was designed 
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as the foundation however, when the structure was redesigned for isolation, W36xl70 
beams were used to connect the columns at the base of the structure. 

Because the isolated and non-isolated structures are two distinct eigenvalue 
problems, due to the free-free nature of the isolated structure, vertical motion, not 
previously seen in a clamped structure, is no longer constrained. The addition of vertical 
displacements at the base is known as uplift. This uplift motion can increase local 
stresses as well as decrease the designed overturning moment. The effects of uplift can 
be seen at some or all frequencies and can have an adverse effect on many types of 
passive isolators. 

B. EARTHQUAKE SAMPLING 

Seismic events occur worldwide and on a frequent basis. Generally, earthquakes 
who’s magnitudes are greater than five can cause significant damage. The following 
uncorrected ground acceleration plotted in Figure 2 is indicative of earthquake 
magnitudes of 6 and above. Various time history data can be obtained through several 
resources such as [Ref. 8]. 
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Figure 2. Imperial Valley Accelerometer Graph From Ref. [8] 

Without formally verifying the enormous earthquake time history data for 
ergodicity, many of the earthquakes observed lasted for less than 30 seconds with the 
majority of the energy focused within the first 10-20 seconds. Time sampling from 
various monitoring stations is generally on the order of lO'^ seconds. Therefore time data 
histories can total to about 6000 data points. To reduce the amount of data and speed up 
the subsequent numerical solutions, correlations can be calculated similar to data signal 
processing. 

The equations presented in this section are outlined in Reference [9] and are 
presented here in summary. This direct characterization is best seen in the spectral 
density function plot. Beginning with the auto-correlation function. Equation (2), of the 
acceleration time history a set of characteristics can be determined from the earthquake 
time history. 


I r 

= lim — J x(t)x(t + T)dt 

(2) 

T 1 0 


Rxx(0) = mean square value 

(3) 

Rxx(ao) = mean value squared 

(4) 
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Although Equations (3) and (4) are useful. Equation (2) is necessary in the 
computation of the power spectral density function. Taking the Fourier transform of 
Equation (2) a one-sided spectral density function is developed, Equation (5). This 
function is better known as the power spectral density function (PDF). 

G^(/) = 2?i?^(T)c-2^Vr (5) 

0 

Plotting Equation (5) versus frequency will yield a graphical representation of the 
energy distribution of the system. There are many algorithms that take advantage of time 
averaging of the sample data to primarily suppress bias errors in the estimate. From the 
PDF an appropriate sampling frequency can be determined. 

To avoid statistical sampling errors the minimum sampling frequency should be 
equivalent to the Nyquist frequency. 

fresolution ( 6 ) 

To increase the resolution and to avoid aliasing the general rule of thumb for the 
sampling frequency is: 

(7) 

Given the value of Equation (7) and the data obtained from plotting Equation (5) 
versus frequency, the previously mentioned 6000 data points can be dramatically 
reduced. Using MATLAB’s power spectral density function on the El Centro earthquake, 
of 1940 time history data reveals it’s distribution plot in Figure 3. From Figure 3 it is 
apparent that the majority of the energy is well below 1 Hz. Therefore, using a sampling 
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frequency of 10 Hz, the time step equivalent of .1 seconds, reduces the proposed 6000 
data points by a factor of 10 to 600. 



Figure 3. El Centro 1940 North South Power Spectral Density Function of the 
Displacement Time History Data 

C. NONLINEAR BASE ISOLATION 

Since many spectral density plots of seismic events indicate a dispersed 
frequency range, viscous dampers aren’t as efficient in energy dissipation than isolators 
that are frequency independent. Frequency independent isolators maintain their unique 
hysteretic loop across a wide range of frequencies. The general topic of base isolation is 
well researched and discussed in [Ref. 10] and is broadly presented in this sub-section. 

The general nonlinear equation that describes general frequency independent 
isolators is: 

F = -Kg(\ + iS){x-Ug) ( 8 ) 
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Where Kg is the dynamic stiffness, 5 is described as the loss factor, x the base motion. 


and Ug is the ground base motion. Graphically, Equation (8) spans the hysteretic loop as 


shown in Figure 4. 


Restoring 

force 



Figure 4. Hysteretic Loop of an Elastic-Pure Plastic Isolator 
The process by which the hysteretic loop is formed is through cyclical 
deformation along a material’s stress strain curve. Deformation below the tensile and 
compressive yield stresses describes typical linear elastic springs. Deformation beyond 
the initial yield stresses results in plastic deformation. In this region the material resists 
additional stress, due to strain hardening, leading to a secondary stiffness rate. In order to 
generate the hysteresis loop the subsequent unloading from the plastic region occurs 
elastically until it plastically deforms again. 

The basic groups of base isolation in use today are linear base isolators and 
nonlinear energy dissipation isolators. Generally, the isolators are not designed to 
plastically deform but to achieve the similar hysteresis loop through hybrid systems. 
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Typical linear base isolators are elastic springs such as rubber or high damping rubber 
(HDR). Typical nonlinear isolators are hysteretic, visco-elastic, friction, or material 
devices. The preferred isolators in use today are the hysteretic nonlinear devices and the 
friction devices. 

1. Rubber Isolators 

Rubber used by it self is categorized as a linear base isolator. Its main advantage 
as an isolator is it’s dynamic shear modulus, G, remains relatively constant over a wide 
range of frequencies and temperatures [Ref. 11]. Additionally, it’s dynamic damping 
coefficient, 5, in Equation (8), is relatively constant at frequencies lower than 50 Hz 
[Ref. 11]. Given that the majority of the seismic energy was concentrated well below 50 
Hz, the use of mbber isolators is desirable. 

A nonlinear isolator is produced when different types of rubbers are used. The 
resulting combined stress strain curve for the composite rubber isolator is tri-linear 
[Ref.ll:p. 199]. At moderate shear strains, such as wind loading and service loading, the 
effective spring stiffness is very low. At larger shear strains, seismic loads, the stiffness 
can increase dramatically until it decreases to the third stiffness rate. However, since the 
deformation is contained within the elastic regions of the isolator there is no hysteretic 
loop generated. 

Additional materials or devices built into a rubber isolator can produce a 
hysteretic effect. The plastic deformation rate would dominate the isolator’s property at 
larger deformations while the elastic stiffness of the mbber would control the subsequent 
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moderate cyclical loading and unloading on the isolator. The effect, during seismic 


loading, is generally modeled as a bilinear model as shown in Figure 5. 



Figure 5. Bilinear Hysteretic Loop 

Where Ki and K 2 are the two distinct stiffness constants and the variable, yp, is 
the yield to post ratio parameter. 

2. Elastomeric Isolators 

Elastomeric bearings are engineered from synthetic rubbers. The properties of the 
bilinear hysteresis loop can be readily designed and manufactured with the use of 
synthetic rubbers. These devices, known as High Damping Elastomeric Bearings 
(HDEB), are frequently used in both structural and industrial use. A comprehensive 
study conducted in [Ref. 12] validates work on performance and efficiency issues 
between polymer type bearings and High Damping Rubber Bearings (HDRB). The only 
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disadvantage of polymer synthetic rubber isolators is that over time the tensile and tear 
strengths become degraded. 

3. Hysteretic Isolators 

Not including the HDEB or HDRB devices, hysteretic bearings are developed 
using metals or composites such as the lead core rubber isolator. By varying materials 
and manufacturing processes a wide variety of hysteretic curves can be engineered. 
Modeling hysteretic systems is mathematically challenging due to the nonlinear 
governing equations. An accepted mathematical model to Figure 3 above, is the Wen 
nonlinear Equation [Ref. 13]. The equation is very adaptive to various hysteretic type 
bearings. Another model, not used in this paper, is known as the Ozdemir element 
[Ref. 14]. This mathematical model is based on similar principles as the Wen element 
however; it incorporates modifications to allow for Shape Memory Alloy (SMA) device 
modeling. 

D. MATHEMATICAL MODELS 

The general formulations of the mathematical models below are developed in 
previous works outlined in the respective subsection’s reference remarks. The 
presentation of the equations represents general summaries of their work. 

1. Bilinear Element 

The bilinear element’s hysteretic loop is shown above, in Figure 5, and is an 
approximation to the dynamic hysteretic loop as shown in Figure 4 above. The hysteretic 
loop is a function of displacement and velocity across the isolator, primary spring 
stiffness, and the secondary spring stiffness. Program code for the hysteretic loop is 
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formulated using basic logical algorithms, and is executed at the given time step interval 
over the time span. 

In developing the algorithm for the hysteresis loop it is evident that the loop is a 
combination of 4 intersecting linear lines. Additionally, the restorative force is a function 
of the spring constant, displacement, velocity, previous deformation state position, and 
the intercept. Ultimately, to identify which part of the hysteresis loop that the current 
deformation data is on, a set of memorized variables is required. The memorized set of 
variables stores which leg the last deformation state was at and the force intercepts of the 
4 lines. 

As an example, an isolator in equilibrium is subjected to a sinusoidal base 
excitation. Assuming that the excitation begins in the positive direction the initial 
velocity and displacement are positive. Given the primary stiffness as Ki, the force is 
calculated from Equation (9). 

F=K^^x + B, 

(9) 

Since the equilibrium position was at zero velocity and zero displacement the value of Bi 
is zero. If the calculated force from Equation (9) is less than a specified positive 
transition force, Fmax> the intercept, B 2 , for the subsequent leg is simultaneously 
calculated. This is conducted by numerically solving for the linear equation with a slope 
of K 2 that intersects the calculated F and Ax position. When the force, calculated in 
Equation (9), is greater than the transition force, F^ax, Equation (10) is used. 

F=K^^x+B^ (10) 
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Continuing along this leg the slope of the next leg, B 3 , is sinaultaneously 
calculated until the velocity becomes negative and the isolator begins to unload. At 
which point Equation (11) is used. 

F=-K^Ax+B^ (11) 

As the deformation progresses along this leg the next leg’s intercept, B 4 is 
calculated as before. Upon surpassing the negative transition force, -Fmax, Equation (12) 
is used. 

F = -K^^x+B^ (12) 

The loop continues until the end of the time history is reached. Additionally, at 
each time step, the intercepts and a location variable are updated and stored. 
Furthermore, reloading and unloading along the same leg is also possible given the data 
stored at each step. 

The two distinct types of bilinear elements used in this paper are referred to as 
“ideal” and “real.” The real isolator models the basic stress-strain curve, where the 
isolator is allowed to strain harden. Therefore, upon successive reload in either the 
tensile or compressive direction the yield force is increased. The ideal isolator is not 
allowed to strain harden and hence has constant positive and negative transition forces. 

The algorithm of the example above is presented in Figure 6 . Incorporating this 
algorithm into the recursive block-by-block formulation is easily conducted because the 
stored data is passed from block to block. However, since the process within a specific 
block is iterative the final block’s data is not updated until the respective block has 
converged. 
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Figure 6. Bilinear Algorithm 
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2. Wen Element 


As stated earlier, the Wen element is a mathematical equation that accurately 
models the hysteretic behavior of nonlinear isolators. The Wen nonlinear equation takes 
the bilinear hysteretic loop from Figure 4 and generates a more realistic hysteresis. The 
Wen element refines the bilinear element’s hysteretic loop through variation of the input 
parameters, as seen in Figure 7. 



Figure 7. Parameters of the Wen's Model From Ref. [14:p. 355] 

The non-dimensional coefficients A, p, and T| control the shape of the 
hysteresis loop. Parameter A influences the initial loading slope and subsequent 
unloading and reloading slopes. The ratio of ^P controls the unloading and reloading 
curvatures. The effect of ^P>1 is equivalent to increasing the (A) parameter and causing 
the loop to curve outward as seen in Figure 7 where ^>1. The parameter r| controls the 
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curves that join the load and unload curves and is best described as the elastic-plastic 
behavior of the material. 

The general formula for calculating the restoring force is: 

K 

F{t)=a-T x(t) + (1 -I- a)F^.z(,t) (13) 

F. 

Where a-r S{t) is the equivalent linear portion and {\+a)FZit) is the 

5 ^ 

r 

nonlinear portion. The symbol a, is the respective yield to post ratio parameter, and z(t) 
is the dimensionless hysteretic displacement that is defined by the nonlinear equation: 



Where, x (t) represents the strain. Rearranging Equation (14) and substituting the 
variable notation of ^ with yin Equation (14) the equivalent Equations (15) and (16) are 
adopted from [Ref. 13:p. 250]. 

z(t) = -;|x(t)|z(0"''lz(t)l- ;0x(?)|z(r)|” +A 's(t) for n = even (15) 

z( 0 = - z{t)" -fix (0k(0r + A X (0 for n = odd (16) 

A more refined formula for the dimensionless hysteretic displacement z(t), from 
[Ref. 15:p. 79], is used in this study. 

z=-7|j|z|zr’'-y^xIzr-t-Ax (17) 

The combination of the nonlinear parameters from Equation (17) yields the 
maximum value of the restoring force. 
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(18) 


To reinforce the effects of each of the respective Wen nonlinear equation 
parameters graphical representations are presented. An arbitrary hysteretic isolator was 
chosen and a sinusoidal displacement was imposed. Figure 8 illustrates the hysteresis 
loop for the isolator and sinusoidal excitation. 



Figure 8. Wen Element Hysteresis Lx)op for Fy-1, A-2, r|-2, P=l, y=l, 8y=l, a=.01 due 

to Unit Sinusoidal Excitation at .5 Hz 
By varying the constant parameters in Equations (14) and (17), the effects, 
delineated by the dashed lines, are presented in Figure 9. 
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Figure 9. Parameter Change Effects On Wen Hysteresis, a) Increase F,. b) Increase 
A, c) Increase (f d) Increase x], e) Increase by. f) Increase (t, g) Increase y 












3. 


Maxwell Element 


The Maxwell element is another type of nonlinear isolator that is constructed 
from a linear spring and viscous damper in series vice in parallel as in a common shock 
absorber. The formulation was developed in [Ref 14] and is included in the code to 
provide an additional nonlinear element. Although viscous dampers and shock absorbers 
are frequency dependant, the overall nonlinearity makes them sufficiently more efficient 
in energy dissipation then a linear spring. By placing the spring and viscous damper in 
series the isolator becomes frequency independent. 

The nonlinear restorative force equation is: 

F{t)=--F{t) + k5{t) (19) 

c 

E. LOAD DEPENDANT ISOLATORS 

Devices that are used in structural applications not only require horizontal 
stiffness but also vertical stiffness as well. The initial loading as well as uplift loads on 
the isolator plays an important role in establishing the horizontal or shear stiffness of the 
bearing. The broad effects of vertical displacements on rubber and elastomeric isolators 
are well documented in [Ref. 9:p. 170], [Ref. 16], and [Ref. 17] and are broadly 
su mm arized in this section. 

For an isotropic material the loading versus the stiffness function is similar to the 
curve outlined in Figure 11. 
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Figure 11. Isolator’s Properties Versus Vertical Displacements 
For small vertical displacements the change in the isolator loads will be small 
enough that the parameters will be dictated by it’s elastic range. Since the vertical spring 
component of the isolator will generally obeys Hooks law. 


<j=Es (20) 

Where a represents the stress, s the strain, and E the elastic modulus. The properties of 
the isolator are then determined by the vertical displacements. 

Generally, when the isolator has been loaded it has stored elastic energy. Rubber 
is particularly strain sensitive which causes the shear modulus, G, to change according to 
the strains imparted due to the connections to the building and ground bearing, 
represented by Figure 12. 
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Figure 12. Shear Stresses on a Load Isolator From [Ref. 9:p.l42] 

The comparison of the total elastic stored energy caused by the shear forces on 
the top and bottom surfaces to the energy obtained from the average shear strain within 
the isolator will influence the overall effective shear modulus. 

rL=-^Jj(^i + ^l.)dydA ( 21 ) 

Where, A is the cross sectional area of the isolator, and t the thickness. Depending on the 
type of isolator, the values of the average shear strain can be calculated and the resultant 
G could then be iterated since G is also a function of Cc (compression strain). 

Expanding the system to more complex systems such as anisotropic materials 
like composites requires a more detailed knowledge of the property versus loading 
curves. Even a simple model can lead to complex numerical formulations if the 
transition was made into the plastic region. This could easily lead to a hysteretic behavior 
in the vertical component due to the continuous loading and reloading in the plastic 
range. 
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m. EQUATIONS OF MOTION AND NONLINEAR ISOLATORS 


Generating the equations of motion of an arbitrary n-degree of freedom structure 
allows the structure to be described both statically and dynamically in a finite 
mathematical manner. By modeling the structure and directly controlling the number of 
degrees of freedom (DOF), the number of differential equations becomes finite and 
manageable. The left hand side of the equations of motion becomes the representation of 
the elastic-linear stmcture while the right hand side captures all the forces generated 
from the linear and nonlinear attributes of the isolators. 

A. GOVERNING EQUATIONS 

Developing a FE model of a structure is covered in various texts [Ref. 18] and 
the basic structural dynamic formulation is covered in [Ref. 19]. The work in this section 
is a basic summary of work developed in these references. 

Taking a continuous model, which has an infinite number of DOF, and 
converting it to a workable yet realistic mathematical model is the goal of FE techniques. 
In order to determine the discretization of the structure a logical determination has to be 
made in respect to retaining a sufficient number of degree of freedoms. The structure can 
then be describe by a continuous mathematical model and is formulated as a partial 
differential equation (PDE). Avoiding the computational effort associated with a PDE, 
the system is transformed into a set of 2"^* order ordinary differential equations (ODE). 
The transformation is conducted by using elements of lumped-parameter models. The 
formulation is covered in many vibration and dynamic analysis textbooks as in [Ref. 6] 
and [Ref. 19] and is summarized below. Additionally, it is computationally useful to 
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transform the coupled ODE into modal space to achieve uncoupled ODEs. Again, these 
methods are well documented and are used in summary in this formulation. 

For an arbitrary base isolated structure the free-free structure is connected to the base 
isolators as shown in Figure 12. 



n‘floor 


' floor 


3*® floor 


2"^ floor 


n 1 floor 


Figure 12. N-Story Base Isolated Structure 
The structure’s nodal degrees of freedom that connect to the isolator are 


contained in the vector (Cset). (Cset) can contain all six degrees of freedom per node but 
generally the lateral and vertical degree of freedoms are chosen for base isolation. 



The governing equations of motion, assuming proportional structural damping, 
[Ref.4:p. 176], for the coupled free-free structure and the isolators, is presented in 
Equation (22): 

+lC]Jx\ +[ji:L{vL=fe,L (22) 

L Jrtxl I Jnil 

Where, (n) represents the total number of degrees of freedom in the system and 
the forcing function, F^Zit), contains the external isolator forces of length (Cset) and 
zeros for non-coupled degrees of freedom. The general solution to the second order 
differential equation is of the form: 

x(t) = xit\+x{t)^ (23) 

The homogenous solution for particular degree of freedom becomes: 

X = + B cos-JI-Co)J) (24) 

The derivation of the particular solution depends upon the arbitrary forcing 
function. For simple forces such as linear springs the solution becomes trivial through 
the manipulation of Equation (22). However, for nonlinear elements the particular 
solution can be significantly complex. 

A method to directly solve Equation (22) is to transform the second order 
differential equation into two coupled single ODEs. Pre-multiplying Equation (22) by 
[M]’* and moving the stiffness and damping components to the right yields: 

jij=-[Mr'[c]|x|-[Mr'[/s:Xx}+[Mr‘fe} (25) 

The two-coupled single ODEs that represent Equation (25) are: 
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X 


(26) 



|J|=-[ji/r' [c]|i| - wr w*)+ wr te, 1 

The solution to Equation (26) can be easily completed using any validated ODE 
solver, however, it can become significantly problematic if the forcing matrix 
WV {FcZ } is nonlinear. 

B. VOLTERRA INTEGRAL FORMULATION 

The solution of Equation (22) can also be obtained using the Volterra integral 
method. The method of solving for the particular solution of Equation (22) is conducted 
using the convolution integral. The derivation of the convolution integral is developed in 
various vibration textbooks such as [Ref. 6], and is formulated as: 

X (r) = £/( r)h(t - T)dT (27) 

Where, h(t-x), is the impulse response function defined as: 

^(0 = ^ = ^I , sin o)„ (28) 

Where ©„ represents the natural fi'equency, and ^ the damping coefficient. In 
order to sum each discretized time band to yield the entire displacement time history the 
system must be physically realizable: 

^(0 = -t)At= lf(T)h(t - t)dT (29) 

Substituting the particular solution into Equation (23), the complete solution to 
the general equation of motion reduces to: 
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(30) 


The resultant equation is the non-homogenous Volterra integral equation as the 
forcing vector is a (nonlinear) function of the response and its derivative. The 

forcing vector /^'“(t) represents the isolator force due to ground excitation on the DOF 
denoted in the subscript. In order to reduce the computational effort for Equation (30), 
Equation (22) is transformed into modal space. The transformation equation is defined 


as: 




(31) 


Where, n is the number of modes, i is the i* DOF, and j is the j*** mode. 

Pre-multiply Equation (24) with the transpose of the mass normalized mode 

T 

shape matrix, [<h] , the modal space equivalent is: 


(32) 


The modal solution of Equation (32) yields the equivalent Volterra integral form 
of Equation (33): 
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Since, the system of equations was transformed into modal space the impulse 
response matrix is uncoupled and is denoted by the (~) symbol. Transforming Equation 
(33) back into normal space by substituting in Equation (31) yields: 



Equation (34) contains (n + m) unknowns and (n) coupled equations, where (n) 
represents the total number of nodes and (m) the number of nodal values contained in 
(Cset). Clearly, we have more unknowns than equations to easily substitute values in and 
achieve a simple solution. However, since there are only (m) linear or nonlinear forces, a 
recursive iterative solution process can be conducted solely on a reduced (m) sized 
system of coupled equations. The set of reduced nodes, (Cset), contains the associated 
connected base nodes denoted on the subscript of i^“''(^). Therefore, the number of 
equations reduces to (m) and the number of unknowns reduces to (2m). This reduction 
can dramatically reduce the solution time in comparison to methods that solve the entire 
(n) governing ordinary differential equations. 


Therefore, Equation (34) reduces to: 



By employing modal participation methods, the computational effort in solving 
Equation (57) can further be reduced. It is shown, [Ref. 4:ch. 6], that for a system 
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transformed into the frequency domain that for a specified nodal maximum displacement 
to a particular (p) mode is as follows: 


[x I 

I /I max 



(36) 


Where, the superscript identifies the mode, the subscript identifies the node, and (n) the 
total number of modes. From the resultant frequency response spectrum, the 
contributions of the modes closest to the mode in question will have the most significant 
impact. Therefore, the transformation matrix [O], is further reduced from a (n x m) 


matrix into one that is (m x r). Where, (r) represents the number of modes retained. 

The time domain normal mode truncation method approaches the number of 
modes retained in a similar manner. For example, using the i* DOF’s modal 
displacement time history solution, defined by Equation (39), it is obvious that the more 
modes retained the accuracy of the solution inereases. 

(0 + (t) + (0 (37) 

However, truncating to a sufficient number of modes with an acceptable error, can 
reduce the computational effort dramatically. 

Equation (35) now reduces to the reduced Volterra integral form: 

(r) 

rfr(38) 

(r) 

In order to reduce Equation (38) to represent a single i* DOF time history, the 
equation is simplified by multiplying the matrices [0][h][ Multiplying out [<I>][h] 



yields: 


31 






n r 




1^1 




0?^22(^-T) - 








^2 

’ + 


hn it-r) 

^2^22(f~-r) ••• 

<I>5 hrr (/- T) 


•• <I>1] 





X,„ 

h 

0 

,<D* hu{t-r) 



1- 





(39) 


Continuing, the time history nodal response for the i* DOF yields: 




rt 

hjj(t- r)if0{r)<ir+... 

J >1 j=i 


r r 


'L^iKtl'jj(.t-T)l^‘’(T)dz (40) 

J 0 

Where, the notational superscript indicates the respective retained modes and the 
subscript identifies with the DOF values, (i), in (Cset). Therefore, the summation 
notation for (j) signifies that the associated DOFs, contained in (Cset), is summed over 
the total number of modes retained from trancation. 


C. RECURSIVE ITERATION FORMULATION 


By employing modal tmncation and retaining the (Cset) sized set of equations, 
the overall computational effort has been dramatically reduced. The solution of the (m)- 
coupled equations is conducted in a recursive iterative manner. From Equation (40), the 
transient response for a particular nodal degree of freedom is dominated by the 

convolution integrals. The isolator forces, within the respective integrals are 

unknown. Therefore, in order to solve for the Xi(t) nodal degree of freedom time history 
an initial guess is conducted. For a simple linear isolator the linear approximation can be 
represented as: 

Fr = -K,M-x,) (41) 
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Where, Kiso represents the linear stiffness, the subscript (i) the nodal DOF, and Xg 
the ground motion. Substituting all the isolator forces into Equation (40) yields: 


x .(0 = 


r . . 

-x^{ry)dT + ... 

J 0 


r • . 

^ (42) 

J 0 

The solution for any one nodal DOF then becomes a function of all the (Cset) 
unknown DOF displacements. To achieve a solution a recursive interative process is 
employed for each nodal DOF. Therefore, Equation (42) is rewritten as: 


X,. 


'(0 = 


r' r 




hjj(t- T)(-K,Jx‘r(I) - x^it))dT+... 


t'^itihjj{t-t)(-K,JxT(r)-x,(mdT (43) 

J 0 J='^ 

The recursive iteration process uses the values to solve for the x'^'"(t) 

values. The iterative process begins by assigning initial arbitrary values to x'^\t) 
through x^'\t). These values are then substituted into Equation (43) and the resultant 
x“'^(7) values are computed. A tolerance check is conducted between x^‘(t) and 
if th® criteria is not met the x^'\t) values are updated to the A:“'"(r) values. 
The process continues until tolerance is met. An algorithm of the iterative recursive 
method is presented in Figure 13. 
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End 

Figure 13. Basic Recursive Iteration Algorithm 

D. LINEAR AND NONLINEAR RECURSIVE BLOCK BY 
BLOCK SYNTHESIS 

The algorithm used to solve the governing Volterra integral can still be 
computationally exhaustive for large, finely discretized time spans and or large (Cset) 
DOF systems. For physically realizable systems the computational effort can be further 
decreased by recursive block-by-block iteration. The convolution integral can be divided 
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into blocks and summed up to arrive at the total solution. The blocks take advantage of 
the linear properties of the convolution integral. 

By the formulation of the convolution integral, Equation (27), the transient 
response spanned over some period of time t, is the summation of the time histories of 
the arbitrary impulse forces acting at specified times, x. 

r)dT + ...^^ f - r)dr (44) 

Therefore, the incorporation of this linearity of the convolution integral into the 
recursive block-by-block method is presented. In general, each block will synthesize 
both a specified transient time history and an associated forcing function. The recursive 
process requires the use of all previous block’s solutions to achieve the total solution. To 
represent the process Equation (44) can be rewritten to represent discrete block times as: 

[f{T)h{t - T)dT = _[/(0...r,)/2(t - T)dt+ T)dT+... 

f^f(t„...t)h(t-T)dT (45) 

Where, the forcing functions /(O.../,) represent the block’s synthesized force 
over the specified block time span. The recursive properties are not outwardly evident in 
the formulation of Equation (45) but is accounted for computationally for each 
successive block. To illustrate the process a system is synthesized using two blocks. 
Equation (45) reduces to: 

- r)dr = - T)dx + I /(tj - x)dx (46) 
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Computationally, Equation (44) is iterated and converged for the first block. The 
next block then proceeds to solve Equation (44) for the second time block. However, the 
first block’s isolator forcing time history solution has to be extended over the second 
block: 

Here, t) represents the impulse response function that begins at time ti 

vice zero. The time history, x(t), for the extended time block f - t)dz is 

the of the next block. Although it can be generally said that any arbitrary 

can be chosen for each successive blocks, a closer approximation to the solution can 
reduce the iterative process. More importantly, the computed extension is added to the 
final recursive iterated solution of the second block. The recursive block-by-block 
method algorithm is presented in Figure 14. 


^f{T)h(t-t)dr= ^ T)dT+ /(0...i,)A(t,...t-r)<fr (47) 
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End 


Figure 14. Recursive Block-by-Block Algorithm 
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For linear isolators such as the one presented in Equation (41) block-by-block 
iterations and block svuninations are straightforward. Each block follows the algorithm 
presented in Figure 16 and Figure 17. However, for nonlinear forces the recursive block- 
by-block synthesis relies more heavily upon the force solutions of the previous blocks. 

A possible model of a nonlinear isolator is presented as: 

= Mx^-x]),(xg-x'),(xg-x.XK,(t),....) (48) 

Where the superscript * symbol indicates the iterated value achieved from the 
algorithm presented in Figure 16. For the nonlinear isolators presented a Runge-Kutta 4, 
(RK4), ODE solver is used. 

Each successive block requires that the solution to the nonlinear equations have 
the correct initial conditions. Therefore, in order to achieve continuity of the isolator’s 
hysteretic response between blocks, the initial condition is interpolated through a direct 
quadratic approximation through the last three force data points of the previous block, as 
shown in Figure 15. 
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Figure 15. Recursive Interpolation 


E. UPLIFT IN RECURSIVE BLOCK-BY-BLOCK SYNTHESIS 


Beginning with Equation (42) the calculation of uplift requires the retention in 
the analysis of the vertical degrees of freedom. Therefore, if the original retained set of 
nodal degrees of freedom were limited to lateral displacements in (Cset), the new (Cset) 
would include the vertical DOFs associated with the Cset nodes. This would increase 
(Cset) by a factor of two. Including bi-directional excitation with uplift would then result 
in tripling the size of (Cset) [Ref 20]. 

For different isolators the uplift properties and their effects on lateral stiffness is 
well documented in [Ref 21], [Ref 16], and [Ref 17]. For simplicity the uplift portion 
of the synthesis was modeled as a linear spring: 

(49) 

Programming interpolating tables can easily be introduced in the developed 
program if known. Therefore, Equation (40) expands to include uplift synthesis: 

x.(t)= Equation(42) + 

y=i 
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Additionally, since many of the parameters are uplift dependant the following 


isolator models are presented: 


r = -Kl(x,-xJ 


(51) 


= f((Xg -X- ),(xg- x‘),(xg- x*),K’(t), T]. (0,y’(0--) 


(52) 


Where, the superscript symbol, *, represents a now synthesized value as depicted 
in the linear approximation of slope (M) below: 

Paramater' (t) = Parameterr^ + M(x "(0) - x~'(0) (53) 

F. SINGLE ISOLATION TRAPEZIODAL EXTRACTION 


Given an arbitrary excitation displacement time history the resulting forcing 
function can be extracted using a trapezoidal integration scheme. The general formula 


taken from any numerous numerical texts [Ref. 22] is presented as: 

^(0 =~^0i + 2^+2^+ 2 / 4 ... + ^^) (54) 

Applying the formulation to the convolution integral yields: 

a:(0) = ^A<A(0)/(0)) (55) 

x(A») = |a((A(0)/(0) + A(AO/(A/)) (56) 

x(«A/) = ■tif/()i(0)/(0) + 2/i(Af)y(A/) + ...(i(kA/)/(mA/)) (57) 


Rearranging Equations (55) through (57) into matrix form jdelds: 
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The impulse response function matrix is a lower triangular matrix that has zeros 
dominating the main diagonal. This in effect negates the weighting factor of the 1/2. The 
integration scheme then turns into a rectangular integration method with an 0 (h) global 
error. To achieve the same 0(h ) approximation of the trapezoidal scheme the weight is 
placed on f(0). Rearranging the matrices in order to solve for f(t) yields; 
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The solution for the extracted forcing function can be accomplished through 
various numerical methods. The trapezoidal integration scheme is a simple way of 
acquiring the synthesized forces given a non-isolated simple structure. Increasing the 
time step would increase accuracy, but can be computationally exhaustive for large time 
histories. 
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V. RESULTS 


To develop MATLAB code for the RBBIEF method workable and especially 
comparative models were used. A sinusoidal wave and the El Centro North South 
Earthquake of 1940 were used for the base excitation. 

Appendix A presents the validation of the trapezoidal method versus NASTRAN 
and using MATLAB’s ode45 function. Using the same model the recursive block-by¬ 
block method is validated against ode45 using a simple spring. The model is then run for 
an ideal bilinear element and compared to the results of the recursive block-by-block 
method 

In order to demonstrate the advantages and disadvantages of the recursive block- 
by-block method 3 models were synthesized using the 1940 El Centro earthquake data. 
The first finite element model is a 4-story, single square bay, that is 25 feet wide and 
having an inter-story height of 17.5 feet. The structure is modeled using a total of 20 
nodes, 120 DOF, four DOF in excitation, and four DOF in uplift. The second finite 
element model is a 4-story, nine square bays that is have a total width of 75 and having 
an inter-story height of 17.5 feet. The building is modeled using a total of 80 nodes, 480 
DOF, with 16 DOF in excitation, and 16 DOF in uplift. The last model is the 30-story 
structure that is 225 feet by 225 feet at the base and an inter-story height of 17.5 feet for 
a total of 525 feet in height. The 30-story structure is modeled using 2,511 nodes, 15,066 
degrees of freedom, 81 DOFs in excitation, and 81 DOFs in uplift. A summary is 
presented in Table 1. 
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1-Bav 4-Storv 

9-Bav 4-Storv 

30-Storv 

Nodes 

20 

80 

2511 

Degrees of Freedom 

120 

480 

15066 

DOF in exictatlon 

4 

16 

81 

DOF in Uplift 

4 

16 

81 

Base Dimension 

25' X 25' 

75' X 75' 

225* X 225' 

Total Height 

70’ 

70' 

525’ 


Table 1. Parameters of Modeled Buildings 


D. PROGRAM VALffiATION 

1. De-coupled Versus Coupled System 

A similar recursive process is validated in [Ref. 24] and is presented in Appendix B to 
validate the trapezoidal method. Using the trapezoidal method, a 9-DOF spring mass 
system, as shown in Figure 16, is used to present the validity of uncoupling a system. 
The equivalent de-coupled system is illustrated in Figure 17. Using a single block in the 
recursive block-by-block method with a linear spring of 10 Ibf^in as the isolator the 
comparison is conducted. Both programs were run using a time step of .01 seconds for a 
ten second duration with a sinusoidal ground excitation of unit amplitude, oscillating at a 
frequency of .5 Hz. 
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Figure 16. 9-DOF Spring Mass System 



Figure 17. Uncoupled Large Mass 9-DOF Spring Mass System 
















Synthesizing a solution using single block in the recursive block-by-block 


program, Figure 18, represents the base displacement time history. 


Nodal response of first node 



Time (s) 

Figure 18. Base Displacement Time History of 9-DOF System Using Recursive 

Block-by-Block Synthesis 

The trapezoidal results plotted along the curve presented in Figure 18, so absolute 
differences Avere calculated between the two data sets and revealed that the maximum 
difference is on the order of 10'*. 
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2. Block-by-Block Synthesis Versus MATLAB's ODE45 


The recursive block-by-block validation process will use the simple 9-DOF 
system is modeled in Figure 16. The process will begin with comparing the two methods 
using a simple linear spring. MATLAB’s ode45 function will use the state space 
representation of the second order differential equation presented by Equation (14). A 
comparison with an ideal bilinear isolator is then made. Because ode45 is an adaptive 
algorithm, memory-type isolators vvith distinct discontinuities at transition points along 
the loop can make the adaptive solver computationally exhaustive. For example, the 
bilinear element has sharp transition points around which the ode45 function will tend to 
oscillate. To avoid these oscillations, additional code was required to force the solution 
to stay within the hysteretic loop. 

Using the same model in Figure 16, an arbitrary (nt) number of blocks were 
chosen producing the same plot as Figure 19. For clarification the number of block (nt) 
was set to three, the linear spring’s value was 10 Ibf/in, and the tolerance in the recursive 
block-by-block synthesis was set at 10’® of displacement and 10"^ for velocity. 
Conducting the ode45 equivalent and plotting the two graphs together depicted no 
appreciable difference. Therefore an absolute difference plot is again conducted in Figure 
19. 
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Figure 19. 9-DOF Spring Mass System, Absolute Base Displacement Differences 
Between a Three Block Recursive Block-by-Block Synthesis and ODE45 
Continuing the comparisons an ideal bilinear element was chosen next. The 
maximum tensile and compressive forces were set to 5 Ibf and -5 Ibf respectively. The 
initial spring constant, Ki, was set to lOlbf with a yield to post ratio of .04. Again, the 
ground excitation was a unit sinusoidal waveform with a frequency of .5 Hz. Both 
programs were conducted with a time step of .01 seconds and a duration of 10 seconds. 
Figure 20, presents the two time displacement histories of the bottom node. The dashed 
line represents the ode45 solution. Figure 21, represents the two hysteretic loops. The 
circles represent the recursive ode45 solution. 
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Ode I Mot 



Figure 20. Three Block Synthesis Versus ODE45 Base Displacements 
Comparisons of a Linear Spring Isolated 9-DOF Spring Mass System, from a Unit 
Sinusoidal Ground Excitation at .5 Hz 



Figure 21. Three Block Synthesis Versus ODE45 Hysteresis Comparisons of a 
Linear Spring Isolated 9-DOF Spring Mass System, From A Unit Sinusoidal Ground 

Excitation At . 5 Hz 
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From Figure 21 it is shown how the ode45 function can oscillate about the 
discontinuous transition comers. Ultimately, Figure 20 indicates a very good 
approximation of the synthesis in comparison to the ode45 solution. 

3. Model Isolator Plots 

A sinusoidal wave of unit amplitude and a frequency of .5 Hz was used on the 
following modeled isolators below in Table 2. 



Properties 

Real Bilinear Element 

K = 500 Ibf/in 

Yield to post ratio = .04 

Max Tensile Force = 50 Ibf 

Max Compressive Force = -50 Ibf 

Ideal Bilinear Element 

K = 500 Ibf/in 

Yield to post ratio = .04 

Max Tensile Force = 50 Ibf 

Max Compressive Force = -50 Ibf 

Wen Element 

Fy = 100 Ibf/in 

Max Elastic Displacement =1 In 

A=8 

n = 2 
g = .0001 
v=.0001 
alpha=.04 

Maxwell Element 

K = 500 Ibf/in 

C = 51 Ibf-s/in 


Table 2. Model Isolator Paramters of Sinusoidal Input 


Figures 22 through 25 plot the respective hysteretic response loops. The duration 
of excitation was ten seconds and at a time step of .01 seconds. The values in Table 2 are 
arbitrary and were chosen to proved illustration in the figures. 
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B. SINGLE BAY 4-STORY BUILDING 


The 1940 El Centro Earthquake Data, Figure 26, displacement time history is 
used for the ground motion for synthesis. The earthquake data was sampled at .1 seconds 
for a duration of 50 seconds. 



Time tsl 


Figure 26. 1940 El Centro North South Ground Motion Time History 

The uplift isolators were modeled as linear springs. The columns and beams are 50 ksi 
steel structural- members with the following specifications: columns, W36x486; first and 
second floor beams, W36xl70; third floor beams are W36xl60; fourth floor beams, 
W36xl50; roof beams, W36xl35. A simple frame illustration of the building is shown in 
Figure 27. The isolator’s properties are listed in Table 3. Figure’s 28-39 
illustrates the hysteresis, horizontal and vertical displacement time responses to the 
earthquake input. The direction of excitation is along the weak axis of bending. 
Additionally, a cutoff frequency of 10 Hz was chosen leading to ten retained modes. 
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Figure 27. Single Bay, 4-Story Building Frame 


II Properties 

Real And Ideal Bilinear 

Elements 

K = 300 Ibf/in 

Yield to post ratio = .04 

Max Tensile Force = 50 Ibf 

Max Compressive Force = -50 Ibf 

K vertical = 30 Ibf/in 

Stiffness Slope = 30 Ibf/in-in 

Wen Element 

Fy = 20 Ibf/in 

Max Elastic Displacement =3 in 

A=3 

n = 4 
g =.0001 
v=.000001 

K vertical = 30 Ibf/in 
alpha=.04 

A,g,n,v,xy,Fy slope = 
.08/in*(A,g,n,v,xy,Fy) 
alpha slope=.008/in*alpha 

Maxwell Element 

K = 300 Ibf/in 

C = 51 Ibf-s/in 

K vertical = 30 Ibf/in 

Stiffness Slope = 30 Ibf/in-in 


Table 3. Single Bay, 4-Story Isolator Properties 





















Isolator Displacement (in) 


Figure 28. Single Bay, 4-Story, Real Bilinear Isolator, Comer Node Hysteresis for the 
1940 El Centro North South Earthquake Ground Motion 


Lateral Comer Node 



Time (s) 


Figure 29. Single Bay, 4-Story, Real Bilinear Isolator, Comer Node Lateral 
Displacement History for the 1940 El Centro North South Earthquake Ground Motion 
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Comer Node Uplift 



Time fs) 


Figxire 30. Single Bay, 4-Story, Real Bilinear Isolator, Comer Node Uplift Displacement 
History for the 1940 El Centro North South Earthquake Ground Motion 
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Figure 31. Single Bay, 4-Story, Ideal Bilinear Isolator, Comer Node Hysteresis for the 
1940 El Centro North South Earthquake Ground Motion 


Lateral Hysteresis of Comer Node 
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Lateral Comer Node 



Time (s) 


Figure 32. Single Bay, 4-Story, Ideal Bilinear Isolator, Comer Node Lateral 
Displacement History for the 1940 El Centro North South Earthquake Ground Motion 



Time Cs') 


Figure 33. Single Bay, 4-Story, Ideal Bilinear Isolator, Comer Node Uplift Displacement 
History for the 1940 El Centro North South Earthquake Ground Motion 
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Figure 34. Single Bay, 4-Story, Wen Isolator, Comer Node Hysteresis for the 1940 El 
Centro North South Earthquake Ground Motion 


Lateral Corner Node 



Figure 35. Single Bay, 4-Story, Wen Isolator, Comer Node Lateral Displacement History 
for the 1940 El Centro North South Earthquake Ground Motion 
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Figure 36. Single Bay, 4-Story, Wen Isolator, Comer Node Uplift Displacement History 
for the 1940 El Centro North South Earthquake Ground Motion 



Isolator Displacement (in) 


Figure 37. Single Bay, 4-Story, Maxwell Isolator, Comer Node Hysteresis for the 1940 
El Centro North South Earthquake Ground Motion 
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Lateral Comer Node 



Figure 38. Single Bay, 4-Stoiy, Maxwell Isolator, Comer Node Lateral Displacement 
History for the 1940 El Centro North South Earthquake Ground Motion 



Figure 39. Single Bay, 4-Story, Maxwell Isolator, Comer Node Uplift Displacement 
History for the 1940 El Centro North South Earthquake Ground Motion 
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Table 4 and Table 5, present the computational results of the synthesis process 


using MATLAB 5.2.1 on a 733 Mhz Pentium III, Dell Computer. 



Blocks 

Element 

4 

8 

16 

24 

Spring 

8.3 

6 

6.4 

8.8 

Real Bilinear Element 

7.6 

5.4 

6.6 

8.6 

Idea! Bilinear Element 

5.5 

4.3 

6 

8.5 

Wen Element 

5.1 

4.6 

6.1 

9 

Maxwell Element 

8.2 

6.4 

7.25 

10 


Table 4. Single Bay, 4-Story, 1940 El Centro Earthquake Recursive Block-By-Block 

Synthesis Time Results in Seconds 



Blocks 

Element 

4 

8 

16 

24 

Spring 

5.50E+08 

9.30E+07 

5.17E+07 

4.22E+07 

Real Bilinear Element 

1.64E+08 

7.11E+07 

4.78E+07 

3.90E+07 

Ideal Bilinear Element 

1.23E+08 

6.00E+07 

4.38E+07 

3.90E+07 

Wen Element 

2.10E+08 

8.89E+07 

4.40E+07 

3.77E+07 

Maxwell Element 

1.50E+08 

7.04E+07 

4.93E+07 

4.14E+07 


Table 5. Single Bay, 4-Story, 1940 El Centro Earthquake Block-By-Block Synthesis 

FLOP Count Results 

From the proceeding tables there is an indication that although the number of 
FLOPS decrease per block the is an optimum number of blocks. This is due to the 
number of optimum iterations per block. Upon converging to a munber of iterations per 
block this number remains relatively the same for each block. Therefore, increasing the 
munber of blocks with the same number of iterations increases the computation time. 
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C. 9-BAY, 4-STORY BUILDING 


Increasing the model’s complexity, a 9-bay 4-story building, Figure 40, is now 
investigated using the 16 hysteretic isolators that include uplift. Uplift is again calculated 
as in the single bay, 4-story structure. 



Figure 40. 9-Bay, 4-Story Building Frame 

The building is design with 50 ksi steel members. The columns are W36x486, 
first and second floor beams are W36xl70, third floor beams are W36xl60, the fourth 
floor beams are W36xl50, and the roof beams are W36xl35. 

















The structure is synthesized for various isolators listed in Table 6, with a base 
excitation of the 1940 El Centro earthquake. The transient time history is developed 
from a time step of. 1 seconds for a duration of 50 seconds. 



Properties 

Real And Ideal Bilinear 

Elements 

K = 300 Ibf/in 

Yield to post ratio = .04 

Max Tensile Force = 50 Ibf 

Max Compressive Force = -50 Ibf 

K vertical = 30 Ibf/in 

Stiffness Slope = 30 !bf/in-in 

Wen Element 

Fy = 20 Ibf/in 

Max Elastic Displacement =3 in 

A=3 

n = 4 
g =.00001 
v=000001 

K vertical = 30 Ibf/in 
alpha=.04 

A,g,n,v,xy,Fy slope = 
.08/in*(A.g,n,v,xy,Fy) 
alpha slope=.008/in*alpha 

Maxwell Element 

K = 300 Ibf/in 

C = 51 Ibf-s/in 

K vertical = 30 Ibf/in 

Stiffness Slope = 30 Ibf/in-in 


Table 6. 9-Bay, 4-Story Isolator Properties 

Figures 41-52, illustrate the results of displacement time histories and hysteretic 
loops of the respective isolators. Again, the cutoff frequency was set to 10 Hz retaining 
20 modes. The direction of excitation is along the weak axis of ha ndin g 
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Figvjre 41. 9-Bay, 4-Stoiy, Real Bilinear Isolator, Comer Node Hysteresis for the 1940 El 
Centro North South Earthquake Ground Motion 


Lateral Corner Node 



Figure 42. 9-Bay, 4-Story, Real Bilinear Isolator, Comer Node Lateral Displacement 
History for the 1940 El Centro North South Earthquake Ground Motion 
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Figure 43. 9-Bay, 4-Story, Real Bilinear Isolator, Comer Node Uplift Displacement 
History for the 1940 El Centro North South Earthquake Ground Motion 



Figure 44. 9-Bay, 4-Story, Ideal Bilinear Isolator, Comer Node Hysteresis for the 1940 
El Centro North South Earthquake Ground Motion 
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Lateral Corner Node 



Figure 45. 9-Bay, 4-Story, Ideal Bilinear Isolator, Comer Node Lateral Displacement 
History for the 1940 El Centro North South Earthquake Ground Motion 



Figure 46. 9-Bay, 4-Stoiy, Ideal Bilinear Isolator, Comer Node Uplift Displacement 
History for the 1940 El Centro North South Earthquake Ground Motion 
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Figure 47. 9-Bay, 4-Story, Wen Isolator, Comer Node Hysteresis for the 1940 El Centro 

North South Earthquake Ground Motion 


Lateral Corner Node 



Figure 48. 9-Bay, 4-Story, Wen Isolator, Comer Node Lateral Displacement History for 
the 1940 El Centro North South Earthquake Ground Motion 
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Corner Node Uplift 



Figvire 49. 9-Bay, 4-Story, Wen Isolator, Comer Node Uplift Displacement History for 
the 1940 El Centro North South Earthquake Ground Motion 



Figure 50. 9-Bay, 4-Story, Maxwell Isolator, Comer Node Hysteresis for the 1940 El 
Centro North South Earthquake Ground Motion 
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Corner Node Uplift 



Figure 51. 9-Bay, 4-Story, Maxwell Isolator, Comer Node Lateral Displacement History 
for the 1940 El Centro North South Earthquake Groimd Motion 


Lateral Corner Node 



Figure 52. 9-Bay5 4-Story5 Maxwell Isolator, Comer Node Uplift Displacement History 
for the 1940 El Centro North South Earthquake Ground Motion 
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Again, Table 7 and Table 8, present the computational results of the synthesis 


process using MATLAB 5.2 on a 733 Mhz Pentium III, Dell Computer. 



Blocks 

Element 

2 

4 

8 

16 

Spring 

124 

59 

50 

71 

Real Bilinear Element 

50 

49 

44 

47 

Ideal Bilinear Element 

61 

38 

39 

65 

Wen Element I 

57 

38 

40 

66 

Maxwell Element 

73 

44 

45 

72 


Table 7. 9-Bay, 4-Story, 1940 El Centro Earthquake Recursive Block-By-Block 

Synthesis Time Results in seconds 



Blocks 

Element 

2 

4 

8 

16 

Spring 

9.27E+09 

2.89E+09 

7.37E+09 

1.54E+09 

Real Bilinear Element 

3.17E+09 

2.05E+09 

7.01 E+08 

145E+09 

Ideal Bilinear Element 

4.06E+09 

1.62E+09 

6.74E+08 

1.45E+09 

Wen Element 

2.70E+09 

1.24E+09 

6.39E+08 

1.22E+09 

Maxwell Element 

5.10E+09 

1.69E+09 

6.93E+08 

1.54E+09 


Table 8. 9-Bay, 4-Story, 1940 El Centro Earthquake Block-By-Block Synthesis FLOP 


Count Results 


As before, there indicates an optimum number of blocks. In comparison to the 
single bay 4-story structure the optimum number of blocks were relatively the same 
although the number of FLOPS increased. 
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D. 30 STORY BUILDING 


The 30-story building, Figure 53, is now investigated using 81 various hysteretic 
isolators that include uplift. Uplift calculations are carried from the previous examples. 
The structure is synthesized for various isolators listed in Table 9, with a base excitation 
of the 1940 El Centro earthquake. 



r Properties 

Real And 
Ideal Bilinear 
Elements 

K = 300 Ibf/in 

Yield to post ratio = .04 

Max Tensile Force = 50 Ibf 

Max Compressive Force = -50 Ibf 

K vertical = 30 Ibf/in 

Stiffness Slope = 30 Ibf/in-ln 

Wen Element 

Fy = 20 Ibf/in 

Max Elastic Displacement =3 in 

A=3 

n = 4 
g =.00001 
v=.000001 

K vertical = 30 Ibf/in 
alpha=.04 

A,g,n,v,xy,Fy slope = 
.08/in*(A,g,n,v,xy,Fy) 
alpha slope=.008/ln*alpha 

Maxwell 

Element 

K = 300 Ibf/in 

C = 51 Ibf-s/in 

K vertical = 30 Ibf/in 

Stiffness Slope = 30 Ibf/in-in 


Table 9. 30-Story Isolator Properties 

The transient time history is developed from a time step of .1 seconds for a 
duration of 50 seconds. Again, excitation is in weak axis of bending and Figures 54-65 
represents the results of the displacement time histories and hysteresis. 
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Figure 53. 30-Story Building Wire Frame 




Figure 54. 30-Story, Real Bilinear Isolator, Comer Node Hysteresis for the 1940 El 
Centro North South Earthquake Ground Motion 


Lateral Corner Node 



Figure 55. 30-Story, Real Bilinear Isolator, Comer Node Lateral Displacement History 
for the 1940 El Centro North South Earthquake Ground Motion 
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Figure 56. 30-Story, Real Bilinear Isolator, Comer Node Uplift Displacement History for 
the 1940 El Centro North South Earthquake Ground Motion 



Figure 57. 30-Story, Ideal Bilinear Isolator, Comer Node Hysteresis for the 1940 El 
Centro North South Earthquake Ground Motion 
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Lateral Corner Node 



Figure 58. 30-Story, Ideal Bilinear Isolator, Comer Node Lateral Displacement History 
for the 1940 El Centro North South Earthquake Ground Motion 



Figure 59. 30-Story, Ideal Bilinear Isolator, Comer Node Uplift Displacement History 
for the 1940 El Centro North South Earthquake Ground Motion 
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Lateral Hysteresis of Corner Node 



.-150 I-1-1-«-«-1-1 

-15 -10 -5 0 5 10 15 

Isolator Displacement (in) 

Figure 60, 30-Story, Wen Isolator, Comer Node Hysteresis for the 1940 El Centro North 

South Earthquake Ground Motion 


Lateral Corner Node 
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.-I I_I_^ _I_I-1-1 
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Time (s) 

Figure 6L 30-Story, Wen Isolator, Comer Node Lateral Displacement History for the 
1940 El Centro North South Earthquake Ground Motion 
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Figure 62. 30-Stoiy, Wen Isolator, Comer Node Uplift Displacement History for the 
1940 El Centro North South Ground Motion 
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Figure 63, 30-Story, Maxwell Isolator, Comer Node Hysteresis for the 1940 El Centro 

North South Earthquake Ground Motion 
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Lateral Corner Node 



Figtire 64. 30-Story, Maxwell Isolator, Comer Node Lateral Displacement History for 
the 1940 El Centro North South Earthquake Ground Motion 


Comer Node Uplift 



Figure 65. 30-Story, Maxwell Isolator, Comer Node Uplift Displacement History for the 


1940 El Centro North South Earthquake Ground Motion 




Again, Table 10 and Table 11, present the computational results of the synthesis 


process using MATLAB 5.2 on a 733 Pentium III, Dell Computer. 



Blocks 

Element 

2 

4 

6 

8 

Spring 

1388 

823 

783 

870 

Real Bilinear Element 

840 

611 

615 

782.7 

Ideal Bilinear Element 

847 

402 

489 

625 

Wen Element 

819 

576 

785 

989 

Maxwell Element 

603 

485 

576 

680 


Table 10. 30-Story, 1940 El Centro Earthquake Recursive Block-By-Block Synthesis 

Time Results in seconds 



Blocks 

Element 

2 

4 

6 

8 

Spring 

1.21E+11 

4.50E+10 

2.91 E+10 

2.34E+10 

Real Bilinear Element 

6.80E+10 

3.23E+10 

2.27E+10 

2.06E+10 

Ideal Bilinear Element 

3.15E+10 

1.97E+10 

1.74E+10 

1.60E+10 

Wen Element 

9.13E+10 

4.50E+10 

2.90E+10 

2.33E+10 

Maxwell Element 

3.82E+10 

2.56E+10 

2.11E+10 

1.82E+10 


Table 11. 30-Story, 1940 El Centro Earthquake Block-By-Block Synthesis Flop Count 

Results 

For the above calculations, the cutoff frequency was set to 1.0 Hz, which retained 
a total of 26 modes. 

Looking at the optimum niunber of blocks versus time. Table 12, tabulates the 
fastest times for the varying numbers of isolators and DOFs. The table also compares 
time synthesis for various computers and processor speeds. The synthesized values are 
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taken from a 50 second time duration of the 1940 El Centro earthquake displacement 
data for the nonlinear Wen element. 




9-Bav 4-Storv 

15-Bav 30-Storv 


8 

32 

162 

Modes retained 

20 

20 

26 

NDOF 

120 

480 

15066 

Optimum Number of 
Blocks 

8 

5 

4 


10 

8 

8 


15 

60 

1145 

Dell Pentium III 750 
MHz Time (s) 

CO 

38 

576 


Table 12. Optimum Times for Nonlinear Wen Element Synthesis 

As Table 12 indicates that for fairly large base isolated structures the optimum 
number of blocks doesn’t deviate much from smaller models. Although relatively 
constant the timesavings can be significant from the original integral method. 
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VI. CONCLUSIONS 


The recursive block-by-block method is considerably efficient in computing the 
seismic transient dynamic response of multiple degrees of freedom structures using 
hysteretic type isolators. As Table 12 outlines, there are limits to the number of blocks 
used to save computation time. Though the number of iterations is generally reduced per 
block, the recursive iterative process still requires a minimum amount of iterations to 
converge. Therefore, by increasing the number of blocks will incur added synthesis time 
even though the number of floating point operations, FLOPS, is decreased. 

Although the initial reductions in time computation are attributed to the retention 
of significant modes and computing the synthesis over the retained (Cset) degrees of 
freedom, it is the recursive block-by-block method that can further reduce it up to an 
additional 60%. For large (Cset) models the major computational elements are the 
impulse response matrix assembly and the recursive iterative synthesis. As stated above 
the iterative synthesis will reach a minimum, however, the impulse response function 
program can contain a significant amoimt of loops which can have a moderate effect on 
time. However, if the time duration and the number of modes retained remain constant a 
flag variable can be entered as to skip this computation if the only thing that changes is 
isolator properties. This can significantly reduce the computation time for large number 
of isolator based structures. 

Since the program only iterates on the retained (Cset) nodal degree of freedoms, 
it will produce the transient dynamic response history and the respective synthesized 
force time history. Because of this, if two structures had an equal number of base 
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isolators and the same number of modes retained, the program’s run time will be similar 
even though structure two had 10 times as many nodes as structure one did. To compute 
non-isolator coupled degrees of freedom time history responses the integral equation 
formulation of the synthesized forces and the associate impulse response function is 
employed. This integral calculation can be computed directly or by the block-by-block 
method. Either method produces the responses quickly and in a matter of seconds. 

The force extraction of an isolated structure using the trapezoidal integration 
scheme proved very accurate. However, the limitations lay in MATLAB’s memory 
intensive backslash command, ‘\’, and the process only works on system that has only 
one node in vibration. 
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VII. RECOMMENDATIONS FOR FUTURE WORK 


The work contained herein was conducted to develop and validate a recursive 
block-by-block program that uses nonlinear hysteretic type isolators for base isolation. 
There have been studies conducted on the exact nature of uplift effects on isolators. The 
shared data from these studies or additional laboratory data can be incorporated into the 
program as data files, which MATLAB can easily interpolate from. Using this data, a 
model can be designed and isolated using a known or fabricated isolator. This model can 
then be used for further validation of the program. Additionally, optimizing the program 
code can greatly reduce the computation time by re-writing the loops into a more 
efficient maimer, specifically, in assembling the impulse response ftmction. Finally, an 
optimization method can be employed to give the designer an optimal isolator for a given 
excitation and set of structural limitations. 
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APPENDIX: TRAPEZOIDAL FORCE EXTRACTION VERSUS 


NASTRAN AND ODE45 

A W36x848 50 ksi steel column was modeled vising 16 nodes spaced equally 
over the total length of the column of 375 feet. The column was chosen from the 
designed 30 story building. The MATLAB code, developed for the trapezoidal force 
extraction’s data was compared to the large mass method used in NASTRAN, [Ref. 23]. 
NASTRAN generated the modal data as well. The excitation was conducted in the weak 
axis of bending of the column. A 9-DOF system was used to validate the extraction 
process against MATLAB’s ode45 flmction. 

The validation is conducted with a sinusoidal input with an amplitude of 1 
centimeter converted to inches with a frequency of .5 Hz. In NASTRAN, the time step 
was conducted at .1 seconds for a duration of 20 seconds. Figure. 66, depicts the 
transient time history of the tip of the beam. Conducting the Trapezoidal force extraction 
program at a time step of .01 seconds and a duration of 20 seconds and a large mass 
representing the ground, yields a solution that is very similar to the one generated by 
NASTRAN. A better representation of the two plots would be to present the relative 
difference, as shown in Figure 67. Because the time response crosses the time axis, 
calculating the relative errors vice differences would generate gross errors. As the plot 
indicates the comparisons are acceptable. 
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Figure 66. NASTRAN Sinusoidal Base Excitation Tip Displacement Results of a 375 

Feet Steel W36x487 Column 









































Using a 9-DOF spring mass system as presented in Figure 16, MATLAB’s ode45 
function was used in accordance with the following state space representation of the 
governing second ordinary differential equation: 

H r-1 

Where: 

/f” = -(io®)Wi)-Abyr] (63) 

in 

Using Figure 17, the trapezoidal extraction method was used using a large mass 
to represent the ground and conducted at time step of .01 for 10 seconds. Figure 68, 
illustrates the time history response of the first mass. Again, the ode45 plot is very 
similar and a difference plot was graphed in Figure 69. From the figures the extraction 
method is an acceptable method in developing a transient response. 
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Figvire 69. 9-DOF Spring Mass System, ODE Versus Trapezoidal Extraction Absolute 

Difference Plot of Base DOF 
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